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The Vector Spherical Harmonic (VSH) model was originally developed to make the 
results of runs of the National Center for Atmospheric Research Thermospheric General 
Circulation Model (NCAR-TGCM) more accessible to scientific investigators. The 
original VSH model was a pure representation of the output fields of particular model 
runs, insofar as no interpolation was performed between model runs. 

At the time investigators had problems accessing the output of the NCAR-TGCM tor 
several reasons. The sheer size of the model meant that the only computers that it could 
run on at the time were supercomputers such as the Cray 1 at NCAR. This meant that a 
special request had to be registered at NCAR to make the model run. But even after the 
model was run. problems occurred because of the size of the output tiles. Depending on 
the length of the run, these could easily reach sizes of hundreds of megabytes, well 
beyond the capabilities of most transfer methods that were available in 1985. Even il 
they could have been transferred, there would have been considerable difficulties in 
storing and accessing these files on the local machines of most investigators. 

The VSH model was first developed to solve these problems. By fitting Associated 
Legendre Polynomials to the output fields of the TGCM the data storage requirement' 
were much reduced. This, in turn, meant that the output fields were much more 
accessible than hitherto. 

The original VSH coefficients were limited to only three fields: the two component' 
of the neutral wind and neutral temperatures. This last field was further limited by 
applying a Bates fit to the vertical. The disadvantage of this approach was that an 
important attribute of the GCM approach to thermospheric modeling was lost: the ability 
to consider cases where diffusive equilibrium was not appropriate. 

Later it was realized that the VSH model was also capable of being a "stand-alone 
model of the thermosphere, provided that suitable interpolations were made between 
model runs of the TGCM and its later incarnation the Thermosphere/Ionosphere General 
Circulation Model (the TIGCM). This model was developed, and the problem of the 
vertical representation was solved by replacing the Bates fit to the temperature with 
vertical splines. A minimal set of basis runs was developed and a version of this "hy brid 

VSH model was distributed to the community. 

This model provided a good representation of the quiet time thermosphere, but it w a' 
very limited in that the only representation of storm-time conditions involved an 



extrapolation of the existing runs to times of high geomagnetic activity. Doing this had a 
number of weaknesses. The first one is that the auroral oval expands during geomagnetic 
storms, and the extrapolation of existing coefficients cannot represent this expansion. 
Second, the thermosphere does not immediately react to geomagnetic storms. There is a 
considerable delay in the density response, in particular, during which the thermospheric 
density increases "ramp up". Similarly, the thermosphere takes time to recover after the 
geomagnetic storm has ended. This report describes techniques that have been developed 
to provide a better thermospheric representation of these storm conditions. But first the 
NCAR-TIGCM and the basic implementation of the VSH model are described. 
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2. The NCAR-TIGCM 


The NCAR-TIGCM (Thermosphere/Ionosphere General Circulation Model) is a 
three-dimensional, time-dependent model of the Earth's neutral upper atmosphere that is 
run on the CRAY-YMP computer at NCAR. The model uses a finite-differencing 
technique to obtain time-dependent solutions for the coupled, nonlinear equations ot 
hydrodynamics, thermodynamics, and continuity of the neutral gas ( Dickinson et al.. 
1981; Roble et al., 1982) and for the coupling between the dynamics and the composition 
( Dickinson et al., 1984). A TIGCM, with a coupled ionosphere and a sell-consistent 
aeronomic scheme, was developed by Roble et al. (1988). More recently, the model has 
been improved by including a self-consistent dynamo model (Richmond et al.. 1992). 
The latest developments in this suite of models include coupling with the stratosphere 
and mesosphere ( Roble and Ridley , 1994). Detailed descriptions of the NCAR models, as 
well as their input parameterizations, have been given in papers by Dickinson et al. 
(1981. 1984), Roble et al. (1982, 1987. 1988), Richmond et al. ( 1992). and Roble and 
Ridley ( 1994). 

Various publications have discussed the use of both the NCAR-TIGCM that is 
described above and the University College London-thermospheric general circulation 
model (UCL-TGCM - Fuller-Rowell and Rees, 1980, 1983; Fuller-Rowell et al.. 1 987 1 
These models have generally proved successful at reproducing the changes that are seen 
in the neutral thermosphere (e.g., Crowley et a!.. 1989; Burns et al.. 1992a; Fuller- 
Rowell et al . , 1994). 

The NCAR-TGCM used a variety of empirical and semi-empirical prescriptions and 
parameterizations. which were discussed in detail by Roble et al. ( 1984) and Roble and 
Ridley (1987). A brief description of these parameterizations is given here. The solar 
heating and photodissociation terms, described by Dickinson et al. (1981. 1984). we re- 
calculated from the Hinteregger (1981) solar EUV fluxes and the Torn et al. < 1980) UY 
fluxes. The Chiu (1975) empirical model of electron density was supplemented by the 
addition of auroral particles using the prescription of Roble and Ridley < 1987). I he 
auroral oval prescription was similar to the statistical patterns desciibed b\ Spito it »/. 

( 1982) and Whalen (1983) and also those given by Feldstetn and Galperin ( 198.''). The 
ion drift velocities were specified by the magnetospheric convection model ot Heelts a 
al. ( 1982). which includes displaced geomagnetic poles (north geomagnetic pole: 78. .A. 
291. 0E; south geomagnetic pole: 74. 5S, 127. OE). 



Although early versions of the VSH model used runs ot the NCAR-TGCM as its 
basis, there were severe limitations in the utility of using this model as the basis of a 
thermospheric density specifications. A large number of inputs to the model used 
empirical paramaterizations rather than first principles calculations, which limited the 
accuracy of the model to the accuracy of these parameterizations. In addition, model 
temperatures were specified as a perturbation about a global mean, rather than as a 
complete calculation. Because the auroral inputs were an “add-on" to this global mean, it 
was very difficult to specify high latitude density and temperature structures adequate!) 
Minor species neutral composition was not included in the TGCM. and it has become 
clear that these gases play a critical role in determining both the thermal and density 
structure of the thermosphere. Similarly, the lack of an interactive ionospheie requirtd a 
reliance on the parameterized Chiu (1975) model in the TGCM. 

In 1988, Roble et al. (1988) developed a version of the TGCM, the TIGCM. that 
addressed these problems. This coupled thermosphere/ionosphere model is a selt- 
consistent coupled Eulerian model of the thermosphere and ionosphere. It utilizes the 
aeronomic scheme developed by Roble et al. (1987) and Roble and Ridley ( 1987). In this 
model the equations describing both the ionosphere and thermosphere are solved on the > 
degree by 5 degree grid mentioned in the previous section. The Heel is et al. ( 1982) model 
is still used to describe high-latitude motion, while the ion drifts at middle- and low - 
latitudes are described by the Richmond et al. (1980) empirical model, as they were for 
the TGCM. Displaced geomagnetic and geographic poles are used with a dipole magnetic 
field. Solar heating is now solved self-consistently once the external source ot EUV and 
UV radiation is prescribed. Semi-diurnal tides have been included as a bottom boundary 
condition. In addition, the minor species (N( : D). N( 4 S), NO. He and An and ion iO + . 
Oi + . NO + , N-> + and N + ) chemistry is calculated and the ion and neutral chemistry 
contributions to the thermal equation are assessed. In addition, the high- latitude ion and 
electron thermodynamic equations are solved tor ion and electron tempeiatuies. 

The primitive equations that describe the NCAR-TIGCM are included in the next 
section. In order to develop the density specification model the TIGCM had to be tested 
and validated as well as developed. Over the years since it was first developed a number 
of studies have been undertaken that have performed this testing and validation. Ciouh ' 
et al. (1989) first compared this model with data measured during solar minimum 
conditions, but most of their comparisons were done with the earlier TGCM. Burraye et 
al. (1992) compared output from the TIGCM with Atmosphere Explorer data at low 
latitudes. Burns et al. (1992a) compared the high-latitude response of a time-dependent 
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run of the model with the equivalent period of Dynamics Explorer 2 data, finding 
reasonable agreement between the two (Figure 2). Bums et al. ( 1992b) looked at both the 
recovery period after a geomagnetic storm and at low latitudes, again finding reasonable 
agreement between the model and data. As part of the Lower Thermosphere Coupling 
Study (LTCS ), Fesen and Roble ( 1991 ) compared data from this campaign with a model 
simulation of the same period. They found that there was fair general agreement between 
the model calculations and observations of tidal structures. The following three diagrams 



Figure 1 :TGCM 
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Figure 2: TIGCM 



Figure 3: TIEGCM 


Richmond et al. (1982) have described this new model, and the following discussion 
is based on their description. In addition to the self-consistent calculations of the dynamo 







electric fields in the TIEGCM, another improvement over the earlier versions of the 
model has also been made. Earlier versions of the model used a dipole magnetic tie Id 
model for the Earth’s magnetic field. The new version of the model now uses a more 
realistic magnetic field model. 

The ionospheric dynamo model uses the calculated neutral wind together with electric 
conductivities derived from the ion and neutral density distributions to compute the 
electric potential at each time step. Some of the essential features of the dynamo model 
are: 


• Geomagnetic field lines are equipotentials, and current may flow between 
hemispheres along these lines at all magnetic latitudes outside of the polar caps (the 
boundaries of which can be varied, but are normally set to 75N and south) 

• A realistic geomagnetic field (the International Geomagnetic Reference Field 1985) 
is used, with calculations being carried out in magnetic apex co-ordinates. 

• The electric potential distribution is externally imposed within each polar cap. and L 
constrained within the auroral regions (set to be between 60 and 75 magnetic latitude in 
each hemisphere) to approach the Heelis et al. ( 1982) model near the polar cap boundary 

• The equatorial electrojet region is externally imposed. 

All of these changes have the potential to create a more realistic ionosphere, which 
may, in turn, also produce a more realistic thermosphere. 

Very large changes in the neutral wind occur in the equatorial region between the two 
versions of the model. In the fully coupled simulation, the upper nightside jet reaches 
wind strengths of 80 m/s as opposed to strengths of 30 m/s in the case where no 
atmospheric dynamo was included. Similarly, in the daytime the upper westward jet 
increased in speed from 30 m/s to 60 m/s. Lower in the thermosphere, the maximum wind 
speed in the nighttime eastward jet increases from 1 10 m/s to 140 m/s. In general, there 
are marked changes between the neutral thermospheric structure when the dynamo is and 
is not included. 
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3. The Equations Used in the NCAR-TIGCM 


The NCAR-TIGCM solves the basic fluid equations in the upper atmosphere. These 
equations include: the continuity equation; the momentum equations; the thermodynamic 
energy equation: the minor species continuity equation: the 0 + transport equation: and the 
electron and ion energy equations. All of these equations are described here briefly. 

The continuity equation for a single thermospheric species may be written in terms ol 
the partial derivative with time of the mass mixing ratio. v Fi > of the i th species. We use 
horizontal spherical coordinates and a log pressure vertical coordinate (X. <t>. z). w here a. = 
longitude, 0 = latitude and z = ln(P 0 /P), P is the pressure, and P 0 <= 50 |iPa) is a 
reference pressure, to write the following vector equation tor the mass mixing ratio ol 
each species ( Dickinson et al.. 1984). 
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The terms on the rights hand side of this equation represent, respectively, the changes in 
the composition due to molecular diffusion, eddy diffusion, horizontal and vertical 
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advection, and chemical production and loss. The I ppj dz tem1 1S mc 'uded m ’ 10 

eddy diffusion calculations to account for variations in the mean molecular mass with 
height which affect the eddy diffusion rate. The vector mass mixing ratio is given by 
= ( 'Pq: ,4*0 ) and % is defined by: 


n m. 



v = the horizontal velocity vector; u = the eastward neutral velocity: v = the northward 
neutral velocity; w = the vertical neutral velocity; T = temperature; n[ = number density 
of the i^ species, and mi is the mass of the i^ species. Other parameters include D. the 
molecular diffusion coefficient given by 


D = 


"p 

00 

m 

p 

i 

H 

0 

t 


(3) 

where D 0 is a characteristic diffusion coefficient at pressure Poo and temperature T no t = 
2 x 10' 1 cm2 s-l), H 0 is the characteristic molecular nitrogen scale height at T 00 = 273 K 
(= 8.63km), J 02 the molecular oxygen photodissociation rate; k is a rate coefficient toi 
three-body (O + O + M — » O 2 + M) recombination of atomic oxygen given by k = 3.S \ 
IO - 30 exp(- 170/T)/T for M = 0 2 ( Johnston , 1968) and k = 4.8 x 10* 33 for M = N': 

( Campbell and Gray, 1973). K(z) is the eddy diffusion coefficient, hi the mean mass. 
P uo is the atmospheric pressure at the ground (= KP Pa), T 0 o = 273K. §ij is the delta 
function, T is the diffusion time scale. The mixing ratio of N 2 is defined as 


h > = 1 - T -4* 

o, o 


(4) 


Photodissociation provides a source of O given by; 


L -J (?) 

and three body recombination provides a source of 0 2 . The alpha coefficients in the tirM 
term on the right hand side of equation 1 are given by: 

“ll =1 (t) 13 + ( 0 I2 -(5) 13) T 2| 

a 22 = " | ^23 + (^21 " ^23)^1 ! 

"12 = ( 0 12 “ ^. 3)^1 

a,3 = (<J>2| - <^23)^2 (6) 


where the subscript 1 denotes 0 2 ; 2 denotes O: and 3 denotes \ 2 and the Ojj aie defined 
by: 
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where Djj is the mutual diffusion coefficients given by Cole grove I 1966). The matrix 
operator L has elements: 
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where Sq is defined by 

£ = 1 m ' 1 dm 
" m irk dz 


( 9 ) 


The L matrix defines diffusive equilibrium solutions through the equation L HP =0. 
Departures from diffusive equilibrium are driven by the hydrodynamic transport and the 
chemical terms in equation 1. 

The next equations of interest are the eastward momentum equation 
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and the northward momentum equation 
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The thermodynamic energy equation is no longer expressed in terms ot perturbations 
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The minor species compositional equation describes the odd nitrogen species (Nth 
N( 4 S) and N( 2 D) 
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The ionosphere is described by the following equations. The 0 + transport equation 
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Where the ion velocity has parallel and perpendicular components 
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The time rate of change of ion and electron temperature is assumed to be small 
compared with the time scale of the TIGCM, so a quasi-thermal equilibiium can hi. 
assumed for these temperatures. That is. the rate at which the ion and election 
temperatures change as a result of other heating processes is sufficiently fast that we can 
ignore heat transport. 

The electron energy equation is given by 
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The cooling rates that occur as a result of coulomb collisions with ions are given by 

-8 ( T e ■ T. ) 
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The electron cooling rates that result from elastic and inelastic collisions with neutrals 
are given by 
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The ion energy equation. 
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These equations provide the basis for the calculation of the thermosphere ionosphere 
system. Numerical methods are used to solve them simultaneously, and the lesulting 
output fields provide the basis for the VSH model. Before describing the numerical 
schemes used in the VSH model, it is necessary to discuss some problems which 
eventually prevented us from using the electrodynamics version of the TGCM. 



4. Problems with using the NCAR-TIEGCM 


It was originally intended to develop the VSH model further using output tiles 
senerated by the NCAR-TIEGCM. Unfortunately, model validation of this version of the 
model did not proceed as rapidly as was anticipated at the time of the original proposal. 

In the first year and a half of the current contract coefficient, files were developed to 
replace the original TIGCM files. All of the requisite quiet-time model runs, except three 
were made. The last two were indicative of one of the problems that still occurred in the 
TIGCM. These three runs were made during solar maximum conditions, and were for 
seomagnetically active times, with cross cap potentials of 90. In these conditions the 
model blew up relatively soon after the beginning of the day. In the past this has indicated 
problems with shear in the convection pattern. The problem gets worse as attempts are 

made to run the model for storm conditions. 

The other problem that was found during this time is more fundamental, and has since 
been solved, but it indicates that far more intensive testing of the NCAR-TIEGCM needs 
to be undertaken before it is used for anything other than a pure research tool. It was 
found that the electron densities at high latitudes in the TIEGCM (but not in the TIGCM) 
were far higher than could be reasonably expected. This has ramifications not only !<>i 
ionospheric specification, but also for the specification of the neutral thermosphere. Ii 
electron densities are two large at high latitudes, two things happen that affect the neutial 
thermosphere. First, ion-neutral coupling rates will increase rapidly altering the advectne 
distribution of the neutral species and through the pressure gradient force the rate .'I 
vertical transport of constituents and global heating rates through adiabatic expansion and 
compression. The second consequence is that the distribution of Joule heating will 
change. Overall specification of the neutral thermosphere will be in error. 

As was stated this particular problem has been solved, but it is indicative that to this 
time adequate testing of the NCAR-TIEGCM has yet to be canied out. Given thcsi. 
circumstances, it was felt that the TIEGCM was not yet ready to be used as an operational 
model that could provide input to the VSH model. 
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5. VSH Formulation 


The VSH model acts on the output fields from the NCAR-TIGCM. These fields are 
fitted spectrally by a set of Associated Legendre polynomials, and sets of coefficient tiles 
are derived for a variety of different geophysical conditions. The output from the model is 
then obtained by synthesizing these different coefficient files to produce output results 
that are appropriate to the conditions that were put into the model. 

The spectral fit consists of using smooth curves to represent the variation of an> 
output variable (such as mass mixing ratio or temperature) with respect to any input 
variable (such as latitude or ap). Each input/output relation is defined by 1) identifying a 
familiar, easily computable, family of curves (the basis functions) and 2) determining the 
best fitting linear combination of the basis functions to represent the actual atmospheric 
variations. Familiar examples of basis functions include the cosine and sine functions in a 
Fourier fit, and the polynomial functions in a least-squares fit. The weights that make up 
the linear combination of the fit - the spectral coefficients - are what is actually stored in 
the coefficient libraries. 

This process of fitting spectral coefficients on the Cray, then reconstituting them in 
VSH, has two important features associated with it. First, the series generating the 
spectral coefficients are truncated at levels that reduce the computer storage requirement' 
but retain the important variations of the output fields. In this way. computers of modest 
storage capabilities can take advantage of the results of a general circulation model over a 
wide range of geophysical conditions. 

The second feature of this process is that a model based on discrete grid point', 
discrete times, and discrete geophysical conditions is converted into a model that i' 
continuous over these same variables. The nature of the interpolation over each of it' 
continuous domains is described in the following sections. It should be pointed out that 
some of the interpolations (such as the variation over solar activity) aie actually cuuicd 
out across several TIGCM runs (rather than within a run). 
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Horizontal 


On a sphere, the spherical harmonics are the most natural basis functions to describe 
the latitudinal and longitudinal variations of a scalar variable. The spherical harmonica- 
are represented mathematically as: 

Y™(0,X) = P™(cos0) e ,mX , IS) 


where the functions P are the Legendre polynomials, q is colatitude, and 1 is longitude. 
The value of a scalar quantity (e.g. mass density or temperature) at a particular location is 
calculated in VSH from the sum: 


q(9.X) = £ a m „ Y>A) 

m.n 


( 19) 


The coefficients a m , n are complex quantities because of the multiplication by e iml . That 
is. the real part of a m .n will be multiplied by cos ml and the imaginary part will be 

multiplied by sin ml. 

Despite the complicated mathematical form of the spherical harmonics, these 
functions are well understood and their values are easily calculated. The need tor this 
level of sophistication arises because of the convergence of the longitudinal lines at the 
poles (where singularities can occur). The use of spherical harmonics allow the output 
fields to retain continuity across the poles. 

For vector variables, a further consideration is required. The coordinate system is 
singular at the poles, because north and south wind directions are not continuous at the 
poles. A change of coordinates is carried out by expressing the wind as the sum of a 
sradient with zero curl (the velocity potential x) and a curl with zero giadient it Ik 
streamfunction y). That is. the horizontal w'ind field can be gi\en as 

U(0,\) = (0.X) + V x k V|/ (0.A.) ( 20) 

In spherical coordinates, this vector equation is equivalent to the set 
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where r is the radius of the earth. 

The velocity potential x and the streamfunction v are scalar variables that can he 
represented as sums of spherical harmonics using ( 19) and (20): 

;(0.A) = r £ b m . n Y”(0,X) 

m.n 

( 22 ) 

m ^ 

= r 2* c m.n Y n 

m.n 

(The radius of the earth r is introduced as a normalization factor.) Using the identities 

3Y _ 9P e imA. 3Y _ imPe" 11 ^ (23) 

30 30 3X. 

and inserting (23) and (24) into (22): yields the formula for recreating winds trom vector 
spherical harmonic coefficients bm.n , c m.n: 



The distinction between scalar variables (such as temperature or density) and vector 
variables (such as neutral wind) is summarized in the table below. 


Variable 

Scalais 

Vectois 


Basis 

p n imX 
r n c 

m P gimx 3P e im>. 
sin 0 30 


That is, scalars use a scalar basis function, while vectors use a vector basis. The function' 


r> m ^D m 

aEn-. andEU 


are precomputed. 

The coefficient values, b and c, that go into the libraries were computed by performing 
a least-squares fit to minimize the squared error in the spectral representation. TIk 
algorithm for performing the least-squares fit is described in Sdmwztrauber ( 1981). The 
actual values of the streamfunction and velocity potential are normalized by a division by 
V n( n+ 1 ) where n is the meridional index. 
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Vertical 


The vertical structure of a geophysical field is described by the least squares tit ol a 
spline (a piecewise polynomial). A spline of degree n is defined as a piecewise 
polynomial of degree n. for which each of the first n-1 derivatives is continuous. 
Therefore, a cubic spline is a function that is continuous and for which its first two 
derivatives are continuous. The coefficients that define the cubic will generally vary from 
one subinterval to another. 

The motivation for the use of splines to represent vertical variations, is that different 
physical processes occur at different altitudes. As a result, the vertical profiles \ar\ 
greatly with altitude. 

In the TIGCM. the vertical variable is log pressure rather than geometric height 
Therefore, the vertical representation of each of the dependent variables is carried out as a 
function of log pressure. Geometric height is also a dependent variable (as a function ol 
log pressure) with its own set of coefficients in the coefficient library. Therefore, w hen a 
VSH user inputs an altitude, the spline defining altitude versus log pressure is inverted to 
convert the altitude to a log pressure. This log pressure can then be used to retrieve any ol 
the desired output fields. The inversion of the altitude spline is carried out via a New ton - 
Raphson method of finding real roots. 

In the VSH model, cubic splines are used in the altitude range of 1 10-500 km at solar 
maximum. Above about 500 km at solar maximum, the MSIS model is called. The 1 Id- 
500 km range is divided into several subintervals, in each of which the polynomial has 
different coefficients. Each end point of the subinterval is known as a knot. The altitude 
range is different during solar minimum conditions due to the different heights of the 
pressure surfaces, but the principles behind the vertical tit are much the same. 


Temporal 

For both the storm and quiet cases, the spherical harmonic lit described above i- 
performed at every hour of model time. A temporal fit then needs to be made. No spectral 
representation is possible in the storm case so all times have to be used. But the quiet 
time TIGCM runs produce diurnally reproducible results. A temporal representation can 
be used to capture the variations of the spherical harmonic coefficients over the course ol 
a day. A Fourier time series is used for this purpose. 
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A Fourier representation implicitly assumes that the output fields are periodic: that is. 
it is assumed that the fields reproduce themselves on a daily basis. This assumption is not 
valid during times of transient magnetic storm activity. The storm subset of VSH does not 
perform a spectral fit over the course of a magnetic storm, but instead records information 
for each hour of storm time. 

The Fourier fit for the daily variation of any field is represented as follows: 




Ul k=l 

Seven coefficients are retained (T=3). This represents a constant term plus three < time * 
symmetric coefficients plus three (time) antisymmetric coefficients (diurnal, semidiurnal, 
and terdiurnal). 


a, cos kt 
kmz 


TIXlZ " 1 



6. Formulation of a “Storm” VSH model. 


In the original VSH formulation, a set of coefficients was developed for quiet 
geomagnetic times. Later it was realized that such a set of coefficients was inaccurate at 
just those times when thermospheric models were most needed to provide specifications 
for various fields: during geomagnetic storms. Therefore it was decided to develop a 
time-dependent version of the VSH model to create a complete VSH model that was 
applicable to a w ide variety of geophysical conditions. 

Instead of using just one value of Ap. the time-dependent VSH model uses a 24 hour 
time history of three-hourly ap values, running from a time 24 hours prior to the time ol 
observation up to the time of observation. Many problems occur in developing a model to 
simulate conditions during geomagnetic storms, not least of which is the problem of the 
time of storm onset. This introduces another degree of freedom into the problem. In the 
old version of the VSH model, time is represented as a set of Fourier coefficients that 
describe a diurnally reproducible run. Now, with the storm version of the model, this 
quiet time UT dependence of the location of the geomagnetic pole not only has to be 
taken into account, but the time after the storm start also has to be considered. One 
possible solution to this problem would be to make a storm-time run tor each one or two 
hours of storm time and use this to drive the model. However, doing this increases the 
storage requirements of the model considerably and still does not deal with some 
outstanding problems, such as the duration of the storm, and the need to normalize the 
TIGCM to produce accurate densities in this case. Rather than this approach, we used the 
results of research done by Prolss and others (e. g.. Prolss. 1980) which showed that the 
response of the neutral thermosphere to geomagnetic storms is much better oriented in 
geomagnetic rather than geographic space. Using this concept, we were able to develop a 
version of the model that used storm-time perturbation additions to baseline. diurnalK 
reproducible cases. 

Before going into the mechanics of the storm-time version of the model further, it is 
worth examining why such a model is necessary, and why the previous technique of 
extrapolating from time-dependent, diurnally-reproducible runs that were developed for 

active conditions were not entirely successful. 

One major feature of geomagnetic storms is that, unlike the magnetosphere and the 
ionospheric convection pattern, the neutral thermosphere takes some time to respond to 
the onset of geomagnetic storms. The time delay for these changes can vary between 



about half an hour and several hours, depending on the prevailing geophysical conditions 
(Ponthieu et al., 1988). Extrapolation techniques do not allow us to incorporate such a 
delay into the model. 

Another feature that occurs during geomagnetic storms is the large transients that arc 
set up early in the storm (e.g.. Mayr et al.. 1984: Burns and Killeen. 1992a). The-w 
transients are produced during storm-time runs of the TIGCM. but. naturally, do not 
occur in diurnallv reproducible runs. Again, this requires a fully time-dependent model. 

The effects of geomagnetic storms on the neutral thermosphere do not terminate 
abruptly. Instead, there is a long “tail-off' as various processes (especially the effects ot 
NO cooling- Maeda et al.. 1989: Burns et al.. 1989) act both directly and indirectly to 
restore the neutral thermosphere to its quiet-time state. A time-dependent model is needed 
to reproduce this “tail-off . 

As well as these physical reasons for improving the ability to model the time- 
dependent response of the thermosphere to geomagnetic storms, there is a programmatic 
one as well. The present models used to simulate thermospheric densities do a 
particularly poor job during geomagnetic storms. Evidence comes tor this in the numbei 
of satellites lost during geomagnetic storms and the breakdown in communications at thN 
time. The latter may be thought to be an ionospheric problem, but. since Matsushita » 
work in 1959. evidence has accumulated that the loss of electrons in the high-latitude I 
region at this time is largely due changes in neutral composition. Thus, it you want to 
develop an adequate model of the ionosphere, you must use a suitable model ot the 
neutral thermosphere that includes all of the composition, circulation and theimal change'- 
that occur during geomagnetic storms. The time-dependent VSH model does this. 

One of the problems that occurs in attempting to model geomagnetic storms is that 
the thermospheric response to them is still one of the most active areas of thermospheric 
research. In attempting to develop a parametrical model of the thermosphere at these 
times a number of questions had to be posed that have still not been leseaiched piopeih 
One such question is the way that the thermosphere changes in response to high-latitude 
forcing and the length of time that this takes. Similarly , the way that the theimospheis 
recovers after the ion forcing and the length of time that this takes is also not well 
understood. Figure 4 and Figure 5 show our first attempts to study this situation. In our 
first attempt to solve some of these problems we assumed that thermospheric 
temperatures, and by implication densities, would be “saturated" in 6 hours (Figure 4>. 
This Figure was produced by averaging all temperature values at latitudes greater than (-><) 
N for winter conditions at solar minimum. Similar studies were done in the Mimmc] 



hemisphere and for solar maximum conditions. In this case, the maximum value of ap 
was 50 and the storm ran for six hours after the onset. The horizontal axis refers to the 
time after the onset of the storm in hours (storm commencement was 10 UT ). The vertical 
axis is the percentage increase in temperature between the storm run and the equivalent 
quiet-time diumally reproducible case. 

As can be seen from this graph, temperatures kept increasing tor an hour or two after 
the storm had ceased. It is plain that temperatures did not “saturate" during a storm ot ft 
hours duration. Therefore, we increased the duration of the storm to 12 hours (Figure 5). 

Now. maximum average temperatures occur some 10 hours after the storm 
commencement, still two hours before ionospheric forcing ended. It is likely that 
“saturation" has occurred and, thus, that this storm duration could be used as the basis tor 
our time-dependent model. 

The recovery time after the storm presents even greater problems. Little work ha 1 - 
been done to categorize this recovery. Hedin et al. ( 1977) looked at a set ot data from the 
AE satellites to categorize many features of geomagnetic storms. They also looked at 
storm-time recovery in a limited w'ay. Their study and a subsequent one by Porter et al. 
(1981). made using the same data, produced a simple empirical relationship that was used 
to formulate the recovery period in the MSIS model. It is likely that the parameterization'' 
used were to simple to model recovery adequately, although this hypothesis has not been 
rigorously tested using independent experimental data. 
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Figure 4. Temperature changes north of 60 degrees for a storm of 6 hours duration. 


This problem was not addressed again with the use of experimental data tor another 
decade. Burns and Killeen (1992 b) used results from an unpublished manuscript as part 
of a review of the thermospheric composition changes that occur during geomagnetic 
storms. They found a relatively “clean” storm, where there were good data from the 
recovery period. By “clean" we mean that the storm was of relatively short duration and 
that the period preceding the storm was cjuiet tor a consideiable time. Thus, 
contamination from previous storms was unlikely to be a problem. Fuithetmoie. the cut- 
off” or cessation of magnetospheric forcing at the end ot the storm was very sharp, a 
result of a strong northward turning in the z component of the interplanetary magnetic 
field at this time. 
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Figure 5. Temperature changes north of 60 degrees for a storm of 12 hours 
duration. 

Perhaps of greatest interest of this study to the development ot a hybrid model ol the 
thermosphere is the duration of the recovery period. It is clear that, within 18 hours ol the 
end of the storm, no decrease in the ratio is present and that 5CFT of the recovery has 
occurred within 6 or 7 hours after the end of the storm. Temperature recovery is even 
faster, suggesting that density recovery will be rapid after the end of the storm. Anothci 
interesting feature of the recovery is that it is not uniform. It occurs first at the poles, and 
then progressively at lower latitudes. 

The TIGCM reproduces the major features of this plot. The recovery time within t Ik- 
model is similar and the tendency to recover first at high latitudes is also modeled 
successfully. Thus, we can be reasonably confident that the time of recovery indicated m 
Figure 5 is similar to that which occurs in the real thermosphere, and that the gio" 
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features of this recovery are well represented. Thus, we were able to use the TIGCM 
recovery scheme in our time-dependent VSH model w'ith little modification. 

Considerations such as those mentioned above determined the structures that we 
would use to develop the time-dependent VSH model. The precalculation ot the 
difference fields that was necessary to produce this model made it easier to apply 
interpolation schemes to the storm effects without modifying the underlying base model 
that was dependent on the other geophysical conditions. 

The basic model structure is that difference coefficient fields are calculated in 
sjeomasnetic coordinates. These difference fields are then used as "add-ons to existing 
coefficient files that supply the underlying, non geomagneticallv forced, background 
thermosphere. There was a minimum requirement of an additional 12 TIGCM runs to 
produce these difference fields, but their modular arrangement makes it easy to omit them 
if time is a more important concern than accuracy. 

The technique of producing these files is outlined in Figure 6. 



Figure 6: The time-dependent VSH Technique. 
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A problem that arose in developing the time-dependent model was how best to 
represent the changes that occur with time. Previously, the VSH model had always used a 
Fourier transform to reduce the number of coefficients needed. The initial version ot the 
time-dependent VSH model also used this scheme. However, it was soon found that this 
resulted in spurious information being included in the model output, so it was decided to 
produce one coefficient set for each hour of storm-time. 



7. Comparisons Between the VSH Models and SETA Data (solar maximum 
conditions). 


The storm VSH model was first compared with SETA (Satellite Electrostatic Triaxial 
Accelerometer) data for a set of 14 randomly selected days from June 1982. On two of 
the days there was almost continuous storm conditions prevailing, while on another day 
there was a long period during which a storm occurred, and a recovery period. 

The results of this study are shown in Figure 7. In this Figure the time-dependent 
VSH model standard deviations are compared with those calculated by MSIS and with 
those calculated by the non time-dependent VSH model. 

Overall performance solar maximum 
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Figure 7. Standard Deviation Comparisons between VSH phase 2 and the time 
dependent VSH model for the randomly selected days mentioned in the text. 
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In this test the time-dependent VSH model did \°Ic better than the non time-dependent 
version of the VSH model. The relatively small difference between the two versions ol 
the model is attributable to rarity of storm events, even during an active period like the 
one modeled here. It will become clearer when we look at the Kp trace that performance 
during storms is superior for the time dependent model. 

We also broke the model down into bins by isolating the model performance using a 
variety of differing geophysical conditions as criteria. Comparisons were then made 
between the time dependent version of the model and the non time-dependent version to 
see where differences occurred. The results of these studies are shown in Figures 8 to I ! 
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Figure 8. Comparisons between the two versions of the VSH model binned In 
altitude. 

In Figure 8. the altitude plot, the greatest improvement occurs at the highest altitude- 
(230-260 km), where there is a 2% increase in density specification accuracy from II' 
to 9%, but there is not much improvement at low altitudes. Presumably, this large 
improvement at high latitudes occurs because storm effects are much more dominant ai 
high latitudes. In addition, apogee for this mission occurred near the south pole, so both 
the high and low altitude bins occur in predominantly polar regions, where the greate-i 
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improvements are expected to lie. The northern bin is at too low an altitude for the 
densities to be affected as strongly by geomagnetic activity as the other bin. 
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Figure 9. Comparisons between the tw'o versions of the VSH model binned by f 10. 


The F10.7 bins also show marked variations in response to the improvements m 
geomagnetic forcing representation that occur with the time-dependent version ot the 
model. However, these variations are an artifact of the small sampling si^e ot the H<> 
binning. Because we have only fourteen days of data, we have only fourteen independent 
F10.7 samples. Both days upon which storms were almost continuous happened during a 
period when the F10.7 values fell within the middle bin. Thus, there was a marked 
improvement in the density representation in the middle bin ( 14-10^ > and little 
improvement in the other two bins. 
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Figure 10. Comparisons between the two versions of the VSH model binned by kp. 

Perhaps the most interesting Figure of these four is Figure 10. which shows the 
improvement of the time-dependent VSH model over the older version at different kp 
levels. Naturally enough, there is little improvement at low kps. as the time-dependent 
model does not affect densities in this kp range. However, the time-dependent model 
should not affect kps in the middle kp range either, yet there is about a 0.5% 
improvement in this range. Presumably, this difference is due to the time-dependent 
model’s calculation of densities for a recovery period after the cessation of geomagnetic 
forcing. The largest improvements occur at high kps. Now the errors in the time 
dependent model’s estimation of density is 9.5%. as opposed to the earlier version of the 
model, which had an error of 14%. This represents a tremendous improvement in densitv 
specification attained by using a true time-dependent model. 
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Figure 11. Comparisons between the two versions of the VSH model binned In 
Latitude. 

The last Figure in this suite (11) shows the latitudinal variation ot densit> 
specification errors. Like the F10.7 and altitude plots, this plot must be treated with a 
little caution. As we stated earlier, satellite apogee occurred near the southern pole, so the 
southern latitude bins represent density changes at the highest altitudes, while northern 
density changes represent those occurring at the lowest altitudes. Therefore, densitv 
specification improvements are greatest in the southern high latitudes, where they are oi 
the order of two percent. 
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. Comparisons between the VSH model and OSS data from AEC (solar minimum 
riditions) 


Comparisons between AE data and the VSH model are more problematical as the data 
id to be pushing the level of reliability that is achieved by the VSH model. The current 
mparisons are those made with the Open Source neutral mass Spectrometer (OSS). 

The data come from all of the long duration storms that occurred during the first year 
AEC operations. They are detailed in the following table. 


dtc 

Duration 

Maximum kp 

vf)R0->7408 1 

21 hours 

6.667 

K)82 

9 hours 

6.0 

1093 

6 hours 

6.0 

1108 

1 2 hours 

6.0 

i l 86->74 187 

42 hours 

8.667 

1204->74205 

30 hours 

7.0 

4286->74291 

138 hours (not continuous) 

7.0 

43 15->743 16 

30 hours 

6.3 


able. Geomagnetic storms that occurred during the first year of AE-C operations 
>r which there was good data coverage. 


Unlike the solar maximum case all of the data come from periods m which 
■omagnetic storms were occurring. Thus, in the following plots, no comparisons have 
;en made with respect to geomagnetic activity as the amount of data available tor such 

imparisons would be slight. 

The first comparison is an overall one, made between the "storm" VSH model and t te 
,d non-storm version of the model. 
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The first such comparison is by local time. Data existed for local times from 0 to IS. 
but not for local times between 18 and 24, so no bar is given for the last local time band. 


Local Time Comparisons 
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Figure 13. Local time Variations at solar minimum 


The largest uncertainty occurs in the afternoon local time band. This problem relates to 
the nature of the disturbance associated with the geomagnetic storm, and the effects are 
also affected by latitude variations. During storms the whole thermosphere heats up. bin 
the region during the daytime is particularly affected by changes in the pressure gradient 
force, and hence by changes in heating by adiabatic expansion and contraction. The 
NCAR-TIGCM is misestimating these changes during solar minimum conditions, 
although it does not do so during solar maximum conditions. 

The next plot shows the latitudinal comparisons. 
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Latitude Variation 
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Latitude 

Figure 14. VSH OSS differences in terms of latitude. 

Surprisingly, in this plot the greatest uncertainty occurs at low latitudes rather than ui 
the high latitudes where it should be expected during geomagnetic storms. Again this is 
an indication that the NCAR-TIGCM is missestimating the effects of adiabatic 
expansion. We will show shortly that the NCAR-TIGCM is in tact underestimating this 
effect. 

The next plot, which shows latitudinal variations, is somewhat less informative. 
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Altitude Errors 
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Figure 15. VSH OSS differences in terms of altitude above the Earth's surface. 

The density variations are expected to be smaller at lower altitudes where storms cause 
smaller perturbations. However, there may be another effect insofar as the variations due 
to latitude and local time may bend back into the altitude variations, due to the relativeb 
small size of the data sample. 

Can we determine whether the VSH model is overestimating or underestimating the 
storm variations? There is an easy way to do this. Instead ot deteimining the \ at Linn 
between the VSH results and the model data, one uses the mean difference. Figure 16 has 
been developed by doing this. 
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Differences VSH - data 
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Figure 16. 

It is clear that the NCAR-TIGCM is considerably underestimating the densitv 
variations that occur during geomagnetic storms, and that this then altects the VSH 
estimations in turn. 

This complication has been corrected in the following manner. Because we are using 
difference fields, a simple correction factor was applied to the coefficients. The solar 
minimum storm coefficients were increased by the appropriate amount for each latitude 
band. These then fold into the equivalent solar maximum storm coefficients by the 
weighting schemes discussed previously. 




11. By Considerations: Should B y Effects be Included 


As part of the contract it was proposed that the importance of B y (The Y component of 
the IMF) effects should be assessed. Two factors would affect their possible inclusion in 
the model. The first is obviously whether they are important or not. The second is 
whether they are sufficiently important to justify inclusion in the model given that there 
would be a large overhead in terms of model run time and storage requirements. In this 
section it is shown that, in fact, changes in B y during storm time are ot little importance 
for all fields apart from the neutral winds, so they have not been included in the VSH 
model. 

The Interplanetary Magnetic Field (IMF) affects the Earth's thermosphere through 
magnetosphere-ionosphere coupling. Charged particles in the solar wind interact with the 
Earth’s magnetosphere forcing electric fields at high latitudes. These electric tields are 
altered by the strength, and particularly by the direction of the Interplanetary Magnetic 
Field. At present the convection pattern for northward IMF is not well enough known to 
include the results for this condition into models with confidence that it is appropriate 
However this is not a large problem for global thermospheric problems, as the effect ot 
the electric fields is minor when B z is northward, so that approximating northward 
conditions with very low polar cap potentials is sufficiently accurate tor all but the 
highest latitudes. In terms of general thermospheric studies the more interesting case 
occurs when B z is southward (or in certain other circumstances the discussion of which 
requires a level of complexity which is not appropriate here). In this case the sign ot the 
Y component of the IMF (B y ) becomes important. 

Before discussing the effects of B y on the neutral winds and ultimately the neutral 
densities further, it is necessary to describe briefly the TIGCM runs that were made to 
produce these winds. Two model runs are described here. In each case (see Figure 17) the 
same hemispheric powers and cross cap potentials were used. 



40 


Time Variation of Auroral Inputs 
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Figure 17 

Hemispheric power was increased from 11 to 100 GW between 9 UT and 10 LT. The 
model was then run for 12 hours with this power. After 22 UT hemispheric power i' 
decreased steadily from 100 GW to 1 1 GW over an hour. In one case this simulation ot a 
large storm was run for By values of -10 nT, in the other case it was run tor B y values ot 

+ 10 nT. 

Figure 18 shows the resultant winds at 18 UT for B y negative (top) and B v positive 
(bottom). These plots are for the southern (summer) hemisphere, and their outer circle 
occurs at 30 degrees, so they include a large swathe of the middle latitudes. 
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Figure 18 

There are clear differences between the two plots. In the main antisunward channel, 
winds during By positive conditions are rotated by 30 degrees relative to those that occur 
in B y negative conditions. The evening convection cell is much better deimed during B v 
negative conditions than during B y positive conditions, but the morning convection cell i' 
somewhat evident in B y positive conditions, but not during B y negative conditions. 



Difference plots of these two conditions are shown in Figure l l >. 
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Figure 19 

The differences are obtained by averaging the B y positive and negative data and then 
subtracting this average from the B y negative case. The maximum wind differences are ol 
the order of 300 or 400 m/s, quite large values compared with the original wind liekis. 
However these large differences are concentrated in a relatively small area ol the high 
latitudes. As will be shown shortly, these large differences in the winds do not translate 
into significant differences in the other scalar fields. 
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Similar plots were worked out for neutral density. This field is closer to the lields lor 
which the VSH model is required. Figure 20 shows the percentage differences that were 
worked out using the relationship: 

Percent = (DensB> + - Dens ave ) / DenSave * 1 
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Figure 20 

Figure 20 gives these percentage differences in the southern (summer) hemisphere. At 
this UT, 8 hours after the start of the storm, maximum differences in density are only 4' , 
in the summer hemisphere, and that difference occurs only in a very small area near 



midnight. On average density differences between the B y positive case and the average 
storm behavior are negligible. 

The differences are somewhat larger in winter, but their maximum value is still 
considerably less than 10% (Figure 21). 
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Figure 21 

Average values at high latitudes are of the order of about 2%. The maximum changes 
occur near the ‘-throat" region, the region in which the neutral winds associated with the 
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convection pattern converge on the dayside of the auroral oval. The asymmetry of the 
pattern (negative on the left hand side, positive on the right hand side) is due to the 
different directions of the the antisunward winds in the two extreme cases of B v and to 
the different shapes of the two convection cells for these extremes. 

Clearly, at this particular storm time B y effects are of minor importance in determining 
density changes. The remaining question is whether these changes can become important 
at any time during a storm. To test this effect all percentage differences between B, 
positive densities and the average density poleward of 60 degrees were averaged tor each 
hour of storm time. The resultant data are plotted here (Figures 22 and 23) in terms oi 
average percentage difference versus time after the storm commenced. 
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Figure 22 


Figure 22 shows the averaged percentage variation of density in the southern i summer 
hemisphere for B y positive conditions during the storm. It can be seen from this diagram 
that the average variation is very small and that it does not inciease gieatly (it at all 
during the storm (the first 13 hours of the study) or the recovery period. If global or 
hemispheric averages were used the improved accuracy ol the model would be negligible 
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In winter the maximum averaged effect reaches 2% by the end ot the storm peiiod 
Globally (or hemispherically) this difference is insignificant compared with the very large 
changes of neutral density (total neutral densities can increase by 10097 over the high 
latitudes during geomagnetic storms and neutral N; densities can change by several 
hundreds of a percent. Prolss. 1980). 

In this section it has been shown that the changes in density caused by changes in the 
Y component of the IMF are minor compared with the changes that are caused by 
geomagnetic storms themselves. Even without considerations ot the loss ot run time in 
the model, it is not a useful exercise to include information about the Y component ot the 
IMF in the VSH model, with the possible exception of the case when the model is to be 
used as a high latitude wind model. Because the model is to be used by NASA primarily 
as a low latitude density model. B v effects have not been included in the model that is 
delivered here. Instructions for running the model are included in the Appendix. 
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PREFACE 


This manual is divided into four sections. Part l provides an overview of the capabilities of VSH 
This part is oriented towards all users of VSH and provides the minimum information needed b> an end 
user. Part II details how a FORTRAN program would be written to access the VSH routines. It include^ 
examples of main programs. Part III is a technical description of the underlying atmospheric physics and 
mathematical equations used in the VSH model. It is most useful for space scientists and others who arc 
interested in the underlying science that w'ent into the development of the VSH model. The tinal section in 
directed towards a systems analyst or systems administrator. It provides information on storage, run-umc 
requirements, and installation. 
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PARTI: USER’S GUIDE 

CHAPTER 1: OVERVIEW 

1 Description 

The VSH (Vector Spherical Harmonic) Model is a computer subroutine that prov ides a description 
of the composition and dynamics of the thermosphere ( Killeen, et al, 1987). It is capable ot integrating real- 
time data with high-accuracv computer simulations to produce a hybrid model ot the thermosphere. 

If the user wishes to run VSH. he or she specifies values for time, location, and the prevailing 
geophysical conditions. With this input, the following atmospheric variables may be obtained: 

COMPOSITION 

Densities: O. O 2 . N 2 , H. electron, OT total mass 
Mixing ratios: O, N 2 . O 2 


WINDS 

Neutral: zonal, meridional, vertical 
Ion: zonal, meridional 

TEMPERATURE 

Neutral, ion 

PRESSURE 


Two types of input are used to drive the model: the coefficient library and real-time data. I m 
coefficient library contains the output of a large number of runs from a general circulation model ot 1 la- 
thermosphere and ionosphere (the NCAR-TIGCM). The real-time data consists of location, time, solar ilu\ 
and magnetic activity indices, as well as any composition, density, and/or wind intormation as obtained m 
satellite measurement. 

The general VSH approach is to use the real-time information as a refinement to the compute: 

simulations. In this sense, the satellite data is optional. 

VSH can also access the MSIS (Hedin, 1987) empirical model of the upper atmosphere. It a V Ml 
user requests an altitude that is out of range of the TIGCM, the desired information is obtained from ilv 
MSIS model. VSH can be called with altitudes up to 1500 km. although it relies upon MSIS to pro\ nv 
density and composition data above about 500 km. Similarly. VSH can be called for altitudes as low as " 
km. however. MSIS is used below' 100 km. and a mix of MSIS and VSH results are used between inn aim 
I 10km. MSIS is also called to obtain the H density at all altitudes, as this variable is not provided m the 
TIGCM. 
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PART II: PROGRAMMER'S REFERENCE 

CHAPTER 2: WRITING A PROGRAM TO CALL VSH 


1 . Input/Output variables 

To call the VSH subroutines, use the following (or similar) statement: 

CALL VSH(OLTPUT.L r SING.GLON.GLAT,ALT.LTJDAY.Ft07A.FI()7.AP.APHIST.STFLAG.BYl\lF 


OUTPUT VARIABLES: 

Name Type Description 


OUTPUT (25) REAL The array 

1 ) 

2 ) 


3) 

4 ) 

5) 

6) 

7) 

8 ) 
9> 
10 ) 
ID 
12 ) 
13) 
20 ) 
21 ) 
22 ) 

23) 

24) 


f output variables 

Zonal neutral wind(U) 

(positive eastward) 

Meridional neutral wind(V) 

(positive northward) 

Neutral temperature) T) 

Zonal ion drift (Ui) 

Meridional ion drift(Vi) 

Pressure ( 5x l()‘ l ° ;i: e\p( -Zi bars) 
Vertical neutral w ind(W) 

Atomic oxygen mass mixing ratio 

Mass density 

Electron density 

0+ density 

Ion temperature 

Molecular oxygen mass mixing ratio 
Atomic oxygen density 
Molecular oxygen density 
Molecular nitrogen density 
Hydrogen density 

Molecular nitrogen mass mixing ratio 


nits 


K 


K 


m 

m 


. m i 


(Variables 14-19. and 25 are reserved for future use.) 



INPUT VARIABLES: 


Name T y p e 


Description 


Units 


USING (25 ) LOGICAL 


GLON 

GLAT 

ALT 

JDAY 

FI 07 A 

FI 07 

AP 

APHISTC8) 


STFLAG 

BYIMF 


REAL 

REAL 

REAL 

INTEGER 

REAL 

REAL 

REAL 

REAL 


LOGICAL 

REAL 


An array of 25 logical variables (.true, or talse. > 
corresponding to each of the above 25 variables. 

Set a value to .false, to eliminate calculation ot that 
variable (to reduce computation lime). Set all 
desired corresponding output variables to .true. 

Geographic longitude (degrees i 

Geographic latitude ( -90 to +90) (degrees) 

Altitude (90- 1500 km) 'km) 

Julian Day (e.g. 324) 

Previous 90-day average F10.7 solar flux <50-350) 

F10.7 value of solar tlux (50-350) 

3-hr ap index (2-80) 

8*3-hr ap index (2-100) from the day prior to the 
UT time needed. This array is entered in the time 
sequence of the ap data 

A tlag to switch the storm time VSH tlag on or oft. 

If it is .TRUE, the storm VSH is called. 

The value of the Y component ot the IMF. This is 
not implemented at the moment and should be set 
equal to zero. 


Alternative Inputs: 

I. PRESSURE 

In lieu of specifying an altitude, a constant pressure surface can be entered. This option is invoked 
by specifying a value of ALT between -5.5 (bottom) and +4.5 (top). OUTPUT! 6) then outputs the altitude 
(km) of the inputed pressure. 

Alternatively, in lieu of specifying an altitude, a log density may he specified. This option i- 
invoked by entering a natural logarithm between -18 (corresponding to a density ot 10' s ) and -41 
(corresponding to a density of l()" ls ). OUTPUTC6) then outputs the altitude (km) ot the inputed density 


2. TIGCM RUN # . , X1 

In lieu of specifying the geophysical parameters JDAY. FI07A. F 1 07. and AP. a spec i he TIC tC M 
run can be selected by setting JDAY= -run#, where run# is a number 1 - 38. When this option is selected, 
the F107A.F107. and AP selections are overridden. For example. JDAY= -2> would select run 
Descriptions of TIGCM runs are given in Chapter 7 of this manual. Note that this option does not work 

correct l v if STFLAG is TRUE. 


Optional Data Statement Inputs: 


Name Type Description 

[RES INTEGER Controls resolution level of model. 

For IRES = 1 model is much taster, but less 
accurate. For IRES = 5. model is slower, but more 
accurate. The detault value lor IRES is 4. treter to 
chapter 8.) 


This can be input by including the common block RESO in the driver program. 

COMMON/RESO/1RES 


I .0 
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2 Linking VSH with a main program 

To use the VSH module a main driver program must be written, compiled and then linked with the 
VSH subroutines and libraries: VSH. STORM. VSHIMF VSHLIB. MS IS. On a VAX. the appropriate 
command tor linking these tiles is: 

LINK PROGRAM VSH, STORM, VSHIMF. VSHLIB.MSIS 

where PROGRAM is the name of your mam program. 


On UNIX, the appropriate sequence of commands is: 

f77 vsh. 
f77 vshlib.f 
f77 storm. f 
f77 vshimf.f 
f77 msis.f 
f77 gws3.f 

f77 <your program>.f vsh.o storm.o vshimf.o vshlib.o msis.o gws3.o 

where <your_program> is the name of your mam program. VSH can then be run by ty ping 

a.out 

where a.out is the final executable image 
3 Sample Main Programs 


Two sample main programs that call VSH are shown below. The tirsi program generates a polat 
wind field and is similar to the SAMPLE program that is included with VSH. The second program 
compares the total mass density w ith that obtained from a satellite pass. 


PROGRAM WINDS 
C Creates a polar wind field 

DIMENSION OUTPUT! 25) 

LOGICAL USING! 25) 

LOGICAL STFLAG 
REAL APHIST!8) 

DATA USING/2*. TRUE.. 23*. FALSE./ 

DATA STFLAG/.TRUE./ 

DATA ALT/300/. UT/5/. FI07A/I40/.FI07/180/. AP/1 1/. JDAY/1S2/. BYIMF'O, 

DO l = 1.8.1 
APHIST(I)=1 1. 

IFfI GT. 6.)APHIST(I) = 50. 

ENDDO 

OPEN (15. NAME= WINDS.DAT. STATUS= 'NEW') 

WRITE ( 1 5.*) LATITUDE LONGITUDE U V " 

DO 200 LAT= 40.50.5 

DO 1 00 LONG= 90. 1 80.30 

CALL VSH (OUTPUT. USING. FLOATi LONG ). FLOATiLAT i. ALT. 

1 UT. JDAY. FI07.A. F107. AP.APH1ST. STFLAG. 

2 BYIMF) 

WRITE 1 15.900) LAT. LONG. OUTPUT! I ). OUTPUTS) 

100 CONTINUE 



200 

CONTINUE 


STOP 

9(X) 

FORMAT (2I9.2F8. 1 ) 


END 


PROGRAM DENSTRAK 


C Calculate density along a satellite orbit 
DIMENSION OUTPUTi 25) 

LOGICAL USING (25) 

LOGICAL STFLAG 
REAL APHISTI8) 

DATA STFLAG/.TRUE./. BYIMF/0./ 

REAL LAT.LONG 

DATA USING /8*.FALSE...TRUE.. 16*. FALSE./ 

OPEN (I. FILE= SAT.DAT’.STATUS= OLD ) 

OPEN (15. FILE= DENS.DAT. STATUS= 'NEW') 

WRITE ( 15.*) ' UT LAT LONG ALT SAT. DENS. VSH DENS’ 

DO I = 1.8.1 
APHIST(I)=I 1. 

IF(I ,GT. 6.)APHIST(I) = 50. 

ENDDO 

DO 1000 1= 1.1000 

READ 1 1 ,5000.END=2000) JDAY.UT. LAT. LONG. ALT. DEN. FI07A. FI 07. AP 
CALL VSH( OUTPUT. USING. LONG. LAT. ALT. UT. JDAY. FI 07 A. FI 07. AP. 
1 APHIST. STFLAG. BYIME) 

WRITE (15.3000) UT. LAT.LONG. ALT. DEN. OUTPUTS ) 

1000 CONTINUE 
2000 CONTINUE 

STOP 

3000 FORMAT (F5. 1 .F6. 1 ,F7 I ,F7. 1 .E 10.3.E 10.3) 

5(XX) FORMAT (I4.F6. 1 ,F7 I ,F7. 1 ,F7. 1.3EI 0.3) 

END 


4 Hints to reduce computational time 


When VSH is used to generate a global or regional field of values (as is the case with the samp' 
program SAMPLE), the latitude loop should always be outside of the longitude loop This will climinat 
the need to recalculate the basis functions (refer to eh. 6). That is. the basis (unctions are mternail 
recomputed each time the latitude has changed from the previous call. It is to your advantage lo lump 
calls with the same latitude in succession. 

For testing purposes, specilv ing a spec i I ic TIGCM run (or the gcophy sical paiametei s ol a spec 1 1 
run) will reduce the number of coefficient files to be read in. and hence also reduce the compulation time 


The parameter IRES w ithin VSH enables users to sellect a faster, although less accural 
of VSH. There are five levels of speed, with I being the fastest, least accurate level, and > 
slowest, most accurate level. The following table gives the average tradeoff between accuracy 
for different values of IRES. Refer to chapter 6 for more details on the accuracy / speed (radeol I 


C. \ C! ' s > i i 1 

hciriL 1 

sliut x|W 



CHAPTER 3: COEFFICIENT LIBRARIES 


There are currently thirteen variables contained in the main coctticiem libraries: 

1 . Neutral wind - divergent component 

2. Neutral wind - rotational component 

3. Neutral temperature 

4. Ion drift - divergent component 

s. Ion drift - rotational component 

6. Height of pressure surface 

7. Vertical velocity 

5. O mixing ratio 

9. Mean molecular weight / temperature 

10. Electron density 

11. 0 + Density 

12. Ion Temperature 

13. CH mixing ratio 

These thirteen fields are used to compute 13 of the variables that are available tor output Th 
remaining 5 output fields are derived fields or are obtained from MSIS. (Currently onlv IS ot the _ 
variables are used.) The same 13 variables are contained in both the storm tstvan and IMF (bwai 
coefficient libraries. 



CHAPTER 4: VSH PROGRAM DESCRIPTION 

1 Program Design 


The VSH model consists of a set of subroutines and does not contain a main FORTRAN program 
(other than the SAMPLE program). Therefore, a mam driving program must he written lor VSH to become 
operational. 

Upon each call to the primary VSH subroutine, the following algorithm is carried out \ la a "crim 
of subroutine calls. The name of each subroutine is listed below followed by its purpose. Each ot the 
indicated subroutines is called from VSH3BY. STORM or VSHIMF. the primary or controlling 
subroutines. 


1. VSHSOLVE 

Solves for a specific field. Calls main computational and ingestion programs 

2. VALIDATE 

Ensures that all needed variables have been assigned values. For example, it geometric height is mpm 
pressure coefficient file is also required. 

3. WHATSNEW 

Determines whether any of the following have changed since the last 
call: latitude, time, altitude, geophysical parameters, required variables. 

4. COEFFSEL 

Determines which coefficient files are to be read and takes weighted averages ot coetlicicnt lilm 
necessary. 

3. TRUNCIN 

Reads in the truncation structure that is used in the variable resolution processing ol \ SH coelticicni- 
6. COEFFIN 

If geophysical parameters or any variables have changed the applicable coelticients are read in. App ■ 
height modifier to altitude versus pressure level to allign TIGCM global average to MSIS global avc;. 
Also modifies TIGCM global average to agree with global averages obtained from the SETA 2 sutclln. 

6. TIMETR 

If coefficients have changed or if the requested time has changed, a Fourier time synthesis m pert. " i 
for the requested time for each variable, altitude, and zonal and meridional wavenumber. 

7. BASIS 

[f latitude has changed, then for the requested latitude. BASIS Ends the values ol the vector -phm 
harmonics P. V. and W for each zonal and meridional wavenumber. 

K. GETPRESS 

Converts altitude or density into log pressure. This is then used in the solution ol the other Melds 


9 . FIND ROOT 

This routine uses a Newton-Raphson method to obtain pressure at a given altitude or density 

10. BVALU2 ifvfl . , 

This is a special version of the spline synthesis routine B VALLE that includes storm and IMF chan 

the altitude of the pressure surfaces. 


ALTTR 



Using the log pressure obtained in ALTTR1. ALTTR2 performs vertical transformations on each \ aiu 
and for each spatial spectral coefficient rn and n. and reduces B-spline coclticients to a single value h 
requested altitude is above the highest pressure surface. ALTTR2 assumes that winds and lempci.i; 
have reached their asymptotic values and assigns the appropriate variables. 

12. SPACTR 

Performs spherical harmonic synthesis at the requested latitude and longitude tor each variable 
winds, it calls vector spherical harmonics V and W and for other variables, calls the scalar "phci 
harmonic P. 

13. DIVTRAN 

Performs the spherical harmonic synthesis for a non-rotational field. 

14. ROTTRAN 

Performs the spherical harmonic synthesis fora non-divergent vector field. 

15. SCATRAN 

Performs the spherical harmonic synthesis for a scalar field. 

16. IONELEC 

Corrects ion drift results. 

17. DENSITY 

Derives densities O. and CK and N mixing ratios. 

IS. RUNPARAM 

This subroutine obtains necessary MSIS data to augment basic VSH model.. 

13. GETMSIS 

Calls MSIS. if necessary, for high altitudes, low altitudes, and H densities. 


Several additional calls arc made in the time-dependent add-on parts ot the model. The> are liMed in ih 
following table. 


Calls in STORM. 


STORM 

A driver routine that calls the subroutines that determine the storm dilterence ti eldv 


2. STMCOEFFSEL 

Selects coefficient files to be used in storm VSH and calculates the weights tor these tiles. 
STMCOEFFIN 

Reads in the storm coefficient files required. 

4. STMTIMESEL 

Inputs correct coefficients for time. 

5. STMBASIS , 

Determines components of the spherical harmonics as a function ot longitude and latitude. 

6. STMALTTR , ... . , 

Solves the spherical harmonics at the specific pressure corresponding to the altitude -clcacd 


STMSPACTR 
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8 . 

9, 

10. 
I 1. 
12 . 
13. 


Calls routines that perform the synthesis of coefficients lor the non-di\ urgent Meld. 
STMDIVTRAN 

Performs spatial synthesis of the coefficients for the non-rotational field. 


STMROTTRAN 

Performs spatial synthesis of the coefficients for the non-divergent Meld. 
STMSCATRAN 

Performs spatial synthesis of the coefficients for the scalar field. 


STMREVTRAN 

Reads in geomagnetic-to-geographic conversion matrix for the scalar Melds. Performs this conversion. 


GTM 

Transformation of geographic location to a geomagnetic location. 


MTG 

Transformation of a geomagnetic location to a geographic one. 


14. TRANSJVEC 

Transforms geomagnetic vectors to geographic ones. 

15. LATMAG 

Transforms geographic location into LT corrected geomagnetic location to acquire the stoim coot hue n 
at the correct place. 

A further set of subroutine calls are mad in the IMF version o! the model. 

1. BYIMFX 

Sets up conditions if the Y component of the IMF is called 

2. BYCOEFFSEL 

Selects coefficient files to be used and calculates weights tor these tiles 


BYCOEFFIN 

Reads in coefficients for the IMF cases 


4. BYCORDSEL 

Matches the UT time with the By case and calls the routines to set up the geomagnetic coordinates 

5. BYTIMESEL 

Transfers from one array to another to make the arrays compatible w uh a later subroutine. 



PART III: TECHNICAL REFERENCE 


CHAPTERS: TIGCM RUNS 


1 Description 


The VSH model is based on the output from a set of runs of the NCAR Thermospheric 
Ionospheric General Circulation Model (TIGCM). This model is run on a Cray supercomputer in Boulder. 
Colorado. The TIGCM solves the three-dimensional Navier-Stokes equations ot fluid dynamics at each 
time step. Grid points are defined every 5 degrees in latitude and longitude and every half scale height in 
altitude. A time step of about 5 minutes is generally used. 

The particular geophysical conditions chosen for the TIGCM runs are representative ol the 
w ide range of conditions that are experienced in the upper atmosphere. As additional runs of the model 
become available, the new data can be easily inserted as new libraries of coefficients in the VSH model. 

Currently, 57 sets of coefficients from the TIGCM are available. There is a basic set ot IS runs 
(#2 1-38) that represent 3 different levels of magnetic activity ( Ap). 3 days within the vear (Julian davsi 
representing seasonal variations, and 2 levels of solar activity (F10.7A and F10.7). The runs are 
summarized below. 

In the table describing RUNS# 21 - 38. the first number in a number pair represents a Rl N# tor 
solar minimum ( low F10.7); the second number is the RUN# for solar maximum thigh F10.7 ). As an 
example, the run corresponding to a low F10.7A and FI 0.7, an AP ot I 1 . and the December solstice is 
RUN #34. 


Ap 


5 


1 l 
32 


RUNS #21-38 


EQUINOX 

JUNE SOLSTICE 

DECEMBER SOLSTICE 

21.24 

27,30 

33.30 

22.25 

28,31 

34.3“’ 

23,26 

29.32 

35.3s 


In addition to the basic 18 runs (#21-38). there is an optional set of 23 additional runs that can Iv 
manually selected by the user. These additional runs are not normally included with the VSH package, but 
are available upon request. RUNS #1-17 represent the annual variation for AP=5 at both solar minimum 
solar maximum. RUNS #39-41 represent seasonal variations at middle values ot the solar cvcle. Rl NS 
#42-44 represent the variation for the solar min-max Fritz peak study. In addition there are X >ets of storm 
quiet difference coefficient files and the same number of IMF difference files. 

In the standard vsh distribution there are 3*13 coefficient files, each representing a different 
variable (refer to chapter 5 ). Each of the main coefficient files contain tits to the basic set ot runs 2 1 
through 38. The other two sets contain the storm and IMF difference files. 
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Density modification 


The mass density outputs from the TIGCM were tound to require lurther adjustments to bring then 
values in line with climatologically averaged satellite data. This step is needed because the TIGCM in quite 
effective in capturing the density variations about a global mean, but less successtul in computing the actum 
mean value. Therefore, for each run. at each pressure level, a mid-latitude comparison was made between 
the TIGCM and MSIS mass density. This difference between MSIS and VSH averages at each ptCNsiiu 
level are used to correct VSH average presure level splines (refer to chapter St. In this wa\ altitude^ me 
mapped to more reasonable pressure levels, enhancing VSH pretormance lor mass densitv calculations 
well as for other output variables). 


The normalization (correction) factors are contained in the file HTMOD.DAT. Each row m tin- 
file contains the correction factors for each TIGCM run. The first column is the run number followed h> me 
corrections to the seven average pressure level splines. 

Similar to the normalization to MSIS values, density values from the SETA 2 satellite lurihui 
normalize TIGCM average pressure versus altitude splines. The normalization (correction) factor- arc 
contained in the files SETAMOD.DAT. Each row in this file contains the correction factors tor cum 
TIGCM run. The first column is the run number followed by the corrections to the seven u\eruge piv'^uu 
level splines. 



67 


CHAPTER 6: INTERPOLATION METHODS 

The VSH model can produce composition and wind data tor a wide range ol geophysical 
conditions and for any location and time. The heart of the VSH model is its library ot spectral coefficients 
The coefficient libraries are obtained by applying spectral fits to the output fields from TIGCM runs. I sine 
these libraries. VSH performs the inverse operation of synthesizing the desired fields from the spectra; 
coefficients. 

The spectral fit consists of the fitting of smooth curves to represent the variation of any output 
variable (such as mass mixing ratio or temperature) with respect to any input variable (such as latitude >" 
ap). Each input/output relation is defined by 1 ) identifying a familiar, easily computable, family ol curve- 
(the basis functions) and 2) determining the best fitting linear combination ol the basis functions !*■ 
represent the actual atmospheric variations. Familiar examples of basis functions include the cosine am: 
sine functions in a Fourier fit, and the polynomial functions in a least-squares tit. The weights that make up 
the linear combination of the fit - the spectral coefficients - are what is actually stored in the coetliciem 
libraries. 

This process of fitting spectral coefficients on the Cray, then reconstituting them in VSH. has two 
important features associated with it. First, the series generating the spectral coefficients are truncated at 
levels that reduce the computer storage requirements but retain the important variations ot the output lieUN 
In this way. computers of modest storage capabilities can take advantage ol the results ot a genciai 

circulation model over a wide range of geophysical conditions. 

The second feature of this process is that a model based on discrete grid points, discrete umc- 
and discrete geophysical conditions is converted into a model that is continuous over these same \ai iablc' 
The nature of the interpolation over each of its continuous domains is described in the lollowing section^ U 
should be pointed out that some of the interpolations (such as the variation over solar activity ) arc actually 
carried out across several TIGCM runs (rather than within a run). 


1 Horizontal 

On a sphere, the spherical harmonics are the most natural basis functions to describe the latiiudiri 
and longitudinal variations of a scalar variable. The spherical harmonies are represented mathematically ^ 

Y„ (6A)= P„ icos0) e ,mX (S.h 

where the functions P are the Legendre polynomials, q is colatitude, and 1 is longitude. The value *-'i 
scalar quantity (e.g. mass density or temperature) at a particular location is calculated m V SH from f 

sum: 

q<e.x,= y a m .„ Y? ( e.D 

m.n (S._i 


The coefficients a m . n are complex quantities because of the multiplication by e im That is. the real pan ■ 
a m n will be multiplied by cos ml and the imaginary part will be multiplied by 'in ml. 

Despite the complicated mathematical form of the spherical harmonics, these t unci tons arc uc 
understood and their values are easily calculated. The need for this level of sophistication arises because ■ 
the convergence of the longitudinal lines at the poles (where singularities can occur). The use ol sphci is 
harmonics’ allow the output" fields to retain continuity across the poles. 

For vector variables, a further consideration is required. The coordinate system is singular at ..: 
poles because north and south wind directions are not continuous at the poles A change ot coordinates 
carried out by expressing the wind as the sum of a gradient with zero curl iihe velocity potential \ i aiui 
curl with /.ero gradient (the streamt unction vl. That is. the horizontal wind Field can he given as 

U(9A) = Vq(8A) + V x ki|/(0./u tsji 

In spherical coordinates, this vector equation is equiv alent to the set 
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u ( 0 , A.) 
v (9, A) 


i + j_ 3y_ 

r sin 9 3A r 39 

_ 1 + 1 

r 30 r sin 0 BX 


(8.4) 


where r is the radius of the earth. 

The velocity potential x and the streamfunetion v arc scalar variables that can be represented as 
sums of spherical harmonics using (8.1) and (8.2): 

4<e,X) = r £ bm.n Y^ie.X) 

m.n 

(8.5) 

v<e.X) = r £ c m ,„ y: (9 X) 

m,n 


(The radius of the earth r is introduced as a normalization factor.) Using the identities 


3Y _ 3P e imX. 3Y _ j m p e 'mA. (8.6) 

39 39 3A 


and inserting (8.5) and (8.6) into (8.4): yields the formula for recreating winds from 
harmonic coefficients h m .n . ^m.n: 


sector spherical 




r 

im P 


" 3P “ 

>v 


iT 



sin 9 


39 





^m.n 

-3P 

+ r 

^ ^ m . n 

im P 


>> e im> ' 

_v_ 

m.n 


L 39 _ 


_sin 9_ 




The distinction between scalar variables (such as temperature or densitv > and 
(such as neutral wind) is summarized in the table below. 


(8.7) 


vector variable'. 


Variable 

Scalais 

Vectois 


Basis 

p n JmX 
r n c 

m P e imX 3P e irT| /c 
sin 9 39 


That is scalars use a scalar basis function, while vectors use a vector basis. The lunuu n 

m -.-m 

„in m P n 


and 3P " 


sin 9 39 

are precomputed. . , 

The coefficient values, b and c. that go into the libraries were computed bv pertorming a U.m 

squares fit to minimize the squared error in the spectral representation. The algorithm lor performing ih. 
least-squares fit is described in Schwarztrauber ( 1 981 ). Th e actual values ol the Mrcamlunct.on am 

velocity potential are normalized by a division by Vmn+1) where n is the meridional index. 
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2 Vertical 

The vertical structure of a geophysical field is described by the least squares tit o! a spline <a 
piecewise polynomial). A spline of degree n is defined as a piecewise polynomial ot degree n. tor which 
each of the first n-1 derivatives is continuous. Therefore, a cubic spline is a function that is continuous and 
for which its first two derivatives are continuous. The coefficients that define the cubic will generally van 
from one subinterval to another. 

The motivation for the use of splines to represent vertical variations, is that different physical 
processes occur at different altitudes. As a result, some variables may have very different vertical profile's 
in the vicinity of 100 km when compared to near 600 km. 

In the TIGCM. the vertical variable is log pressure rather than geometric height. Therefore, the 
vertical representation of each of the dependent variables is carried out as a function ot log pressure 
Geometric height is also a dependent variable (as a function of log pressure) with its own set of coefficient^ 
in the coefficient library. Therefore, when a VSH user inputs an altitude, the spline defining altitude versus 
log pressure is inverted to convert the altitude to a log pressure. This log pressure can then be used to 
retrieve any of the desired output fields. The inversion of the altitude spline is carried out via a Newton* 
Raphson method of finding real roots. 

In VSH, cubic splines are used in the altitude range of I 10-500 km. Above about 500 km. the 
MSIS model is called. The 110-500 km range is divided into several subintervals. The polynomial has 
different coefficients in each of the subintervals. Each end point ot the subinterval is known as a knot. 


3 Temporal 


With the exception of several experimental storm-time runs, the TIGCM produces diurnalh 
reproducible results. At every hour of model time, the spherical harmonic fit described abov e is performed 
The purpose of the temporal representation is to capture the variations ol the spherical harmonic 
coefficients over the course of a day. A Fourier time series is used for this purpose. 

A Fourier representation implicitly assumes that the output fields are periodic; that is. it i* 
assumed that the fields reproduce themselves on a daily basis. This assumption is not valid dining times ni 
transient magnetic storm activity. The storm subset of VSH does not perform a spectral lit over the course 
of a magnetic storm, but instead records information for each hour of storm time. 

The Fourier tit for the daily variation of any field is represented as follows: 


*Vi.r u “ 



T 

I 


a. cos kt 


+ 


T 



<in kt 


( 8 . 8 ) 


Seven coefficients are retained <T=3). This represents a constant term plus three uimei symmetric 
coefficients plus three (time) antisymmetric coefficients (diurnal, semidiurnal, and terdiurnab. 

4. Julian Day 

To account for the variations over the course ot a year, the TIGCM is urn undei identical 
conditions at the two solstices and at an equinox. The solar annual cycle is described by a pure sinusoidal m 
radiation absorbed by the earth. In particular, the variations of any variable over the course ol a year i, 
currently fit using a two-point trigonometric interpolation. That is. a simple sinusoidal curve is assumed lot 
both the annual and semiannual variations. Therefore, trigonometric interpolation between TIGCM runs i, 
carried out. In future versions of the model, a Fourier fit over the annual cycle will be used as the TK i< M 
coefficient library I'iils up. 
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5. Solar Flux 


Both of the solar flux arguments <F10.7 and F10.7a> influence the variation of most variables. 
Fit). 7a (the 90-day average) has the most bearing at lower altitudes. F10.7 expresses the daily transient 
flux, and becomes important at higher altitudes. These two parameters are synthesised into a single value oi 
FI 0.7' using the relation: 

F| 0.7' = F10.7a + k j ( h )*fF 1 0.7-F 1 0.7a) + k 2 (h) * (F10.7 - FI0.7a) 2 


where h is height (or spline knot location), and k|(h) and k 2 (h) are weighting factors obtained from MSIS 
The value of F10.7' is used for all subsequent calculation of geophysical fields. It is always found to he 
between F10.7 and F10.7a. Because the weights vary with height, the ettecitive FI0.7 value will also be 
heisht-dependent. As a result, different linear combinations of TIGCM runs may be used at different height 
levels. As anticipated, the kj factors increase monotonically with height. This reflects the increasing 
importance ot daily solar activity with increasing height. 


The k|(h) and ks(h) factors are contained within the file F107MOD.DAT. Each row in this nk 
contains the linear (k, ) or'quadrat.c fk 2 ) factors for each TIGCM run. The first column is .he run number 
followed by the weighting factors for each height level. 


6 . Magnetic Activity 


The variations over magnetic activity is carried out by a linear interpolation over Ap in the quid 
time case In the storm case, interpolation is made between Storm-quiet difference fields to provide an 
addition to a diurnallv reproducible case (ap is set to I I internally for the diurnally reproducible addition . 
For the scalar fields the difference coefficients are transformed mathematically and then the ditleiencc lid" 
is reconstituted and added to the quiet time output. The winds are treated differently because ol the lack m 
a suitable algorithm for doing a transform on vector coefficients. The wind difference fields me 
reconstituted and then transformed to geographic coordinates. They are then added to the quiet time ouipu: 


7 . Variable Resolution 

A variable resolution scheme is used to improve the computational efficiency ol \ SH The k;c. 
behind variable resolution is to process only those coefficients that will make a non-neglig.ble contr.hu ., . mi 
to the synthesis of the desired fields. In some instances, very tew ot the high ordet Founei and sphuw.n 
harmonic coefficients are required to capture important thermospheric properties. 

VSH provides the user with a choice of five resolution levels. At the highest resolution, the nns.c 
is most accurate but slowest. At lower resolutions, the model is faster, but less accurate. The ;.x atk.biliiy 
five resolution levels enable users to select appropriate spots on the performance-accuracy trade I a.iw 
Generally speaking, the lower the resolution, the greater the smoothing that occurs w uhin the mo. . o . 

'The resolution level is controlled by the FORTRAN variable IRES. The value of IRES has I x 
preset to level 4, corresponding to nearly full resolution. By changing the value ol IRES m the DAI 
statement, the resolution can be changed from I (lowest resolution . to s . highest resolution i. 

For each resolution level, there is i) a set of truncation levels lor the spatial ioi spheiical haimonw 
synthesis: ..) a set of truncation levels for the UT (or Fourier) syntheses, and ... . a set of vhllereiuc 
thresholds. A difference threshold determines whether an input value is significantly vh Mucm u | m la ' 
the previous call. For example, a call to VSH with UT=I2 followed by a call with LT=| Dl wd, . . . 
invoke a new synthesis of the Fourier time representation. For low resolution, the -phci ic.i hamionic ., . 
Fourier syntheses are highly truncated and the difference thresholds are set relatively high. Foi m-J, 
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resolution, the Fourier syntheses extend through much higher values ot m. n. and i in equations iX.Ti and 
(8.8) and the difference thresholds are lower. 

Associated with each of the five resolution levels is an error threshold value. Each error value 
represents the average deviation from the full resolution value. Coefficients that it excluded, would 
contribute an error exceeding the threshold are marked for retention. Any unmarked coefficients are unused 
in the synthesis, producing computational savings. 
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CHAPTER 7: OBJECTIVE ANALYSIS SCHEMES 


l. Density 

VSH allows the user to enter real-time satellite data to lurther improve upon its accuracy r lie 
model has also been tuned to the climatological values of winds and density that have been accumulated 
over many satellite passes. 

Some coefficient files will contain density modifications based on several passes oi i he Sh T\ 
satellite. The slope and intercept used in the vertical representation of density will he moditied so ih.ii 
their simulated global averages correspond to the observed satellite global averages, (also re ter to chap 


2. Winds 

Where available. Dynamics Explorer (DE) wind measurement data have been merged with the 
TIGCM simulation results. The merging process has been accomplished with techniques ot objective 
analysis (Thiebaux and Pedder. 1987) that are commonly applied in meteorological analysis. In umivj 
objective analysis, a value of a field at a given grid point is computed trom a weighted average ot it- 
computed TIGCM value there and observed values in a neighborhood ot that location. The magnitude i»i 
the weight is inversely proportional to the distance of an observation from its site ot prediction. In tho 
manner, the winds contained in several coefficient sets actually consist ot a background TIGCM wind lieu! 
as modified by observations. 
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PART IV: SYSTEMS ANALYST'S REFERENCE 


CHAPTERS: INSTALLING VSH 

I Program and library information 

The VSH model consists of two components: a set of FORTRAN subroutines and a library oi 
numerical coefficients. The model is provided on a standard magnetic tape. It has been saved using the 
VAX/VMS BACKUP command with the volume identifier VSH and the save set name ot YSH.BCK 

The FORTRAN subroutines are VSH. STORM. VSHIMF. VSHLIB, MSIS and SAMPLE. VSH 
contains the primary subroutines of the package, w hile VSHLIB contains the mathematical libraries. MSIS 
is a stand-alone atmospheric model that is called when the altitude is out ot range ot the values contained 
within the coefficient libraries. A sample main program SAMPLE is also included as a demo. 

The numerical coefficient libraries are available in two versions. The V AX/VMS version contains 
the library coefficients already in the necessary binary unformatted torm. In the non-VAX version, the 
coefficients are stored in ASCII unformatted form. Each coefficient tile is named VARXX.WY. whcic 
XX is a two-digit integer and YYY is either ASC (for ASCII files) or BIN (for binary files). The non-VAX 

version is further described in the following section. 

The VAX version is immediately usable, once the files are transferred trom tape to your computet 


2 Use on a non-VAX machine 


The non-VAX version of VSH contains the coefficient library stored in ASCII tormat. A generic 
conversion program is included with this version to convert the ASCII files into binary. This program is 
named VAR_ASCTOBIN. Running VAR_ASCTOBlN will convert each V ARXX.ASC tile to .1 
V ARXX BIN 

The FORTRAN subroutines VSH FOR. VSHLIB. FOR. MSIS. FOR. and SAMPLE. FOR should 
be recompiled on the machine to be used. In addition. non-VAX users may need to change the name ol two 
functions that are intrinsic to the VAX. These two functions convert characters to integer, 1LHAR1 and 
real variables to integers (NINT). These two functions may appear in VSH. VSHLIB. ASC TOBIN and 


BINTOASC. 



CHAPTER 9: STORAGE REQUIREMENTS 

1 Description 


The storajze requirement on disk is essentially determined by the size ot the eoettieient tiles. These tiles .n,. 
about 360K- TOOK each in binary form for the basic runs. If you requested the files in ASCII torm i tor use on a non 
VAX machine), the ASCII files will be found to be about three times larger than these values 
A breakdown of the storage requirements is the following: 
vsh3bv.f 79314 kB 

storm.V 49198 kB 

vshimf.f 23550 kB 


vshlib.f 

msis.f 

2ws3.f 


36035 kB 
50108 kB 
19703 kB 


and the coefficient files: 

varOl ->var!3 5888304 kB 

stvarOl ->stvarl3 8959520 kB 
bvvarOl ->bvvarl 3 672256 kB 


These values are for the sun UNIX system. 



REFERENCES 


Dickinson. R. E.. E. C. Ridley, and R. G. Roble ( 1981 ) A three-dimensional general circulation model ot 
the thermosphere, J. Geophys. Res.. 86. 1499. 

Hedin, A. E. 1 1987) MSIS-86 Thermospheric Model. J. Geophys. Res.. 92. 4649. 

Killeen, T. L., R. G. Roble. and N. W. Spencer ( 1987) A computer model of global thermospheric wnuK 
and temperatures. Adv. Space Res.. 7. 207-215. 

Swarztrauber R.N. (1981) The approximation of vector functions and their derivatives on the sphere. 
S1AMJ. Numer. Anal., 18. 934-949. 

Thiebaux. H.J. and M.A. Pedder ( 1987) Spatial Objective Analysis. Academic Press, London. 295 pp. 



76 


Date 
Mar. 1991 

Mar. 1991 
Mar. 1991 
Aug. 1991 
Aug 1993 


CHANGES FROM VERSION 1 

Change 


Four new variables have been added: electron density. 0 + number density . ion 

temperature, and 02 mixing ratio. 0 + number density and ion temperature are given m 
the coefficient files as the natural logarithm of their actual values. 

The OUTPUT variable array has been increased to 25 from 14. to incorporate the new 
variables and to allow for future expansion. 

The range of inputs for ap has been extended to (2-60). The range ot inputs tor F- in 
has been extended to (50-350). 

FI 0.7a added as an input variable 


The model extension to 90 km has been validated. The storm version ot YSH has been 
included and validated for neutral densities and winds. Extra difference variable tiles arc 
added. The IMF version is included in the package, but is switched ott because no liclJ 
has been validated. Extra difference fields are added. 
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14. Appendix 2. 

Task Report on tasks 1, 2, and 3. 

Nine tasks were listed in the three year period of this grant. It is the first three ot these 
tasks that are of relevance here. These tasks are: 

1) Update the current quiet-time coefficient files to include representations of nun 
made using NCAR-TIEGCM (National Center for Atmospheric Research - 
Thermosphere-Ionosphere-Electrodynamic General Circulation Model) rather than the 
NCAR-TIGCM (Thermosphere/Ionosphere GCM ) runs that had previously been used. 

2) The temperature representation in the VSH model will be revised to allow the 
vertical temperature structure to be determined without the use of Bates profiles. 

3) Comparisons will be made between the present storm-time version of the model and 
data from storm-time orbits of the DE 2 (Dynamics Explorer 2) and AE (Atmosphere 
Explorer - C, D, and E) satellites. 


Tasks 1. 2 and 3 were completed. The NCAR-TIEGCM was run for all of the 
appropriate conditions. However in doing this several problems were encountered duting 
solar maximum conditions. It was found that the model was highly sensitive in these 
conditions and broke down. Upon creating coefficient files for these runs, it was found 
that the performance of the NCAR-TIEGCM was at this stage worse than that of the 
NCAR-TIGCM. so in consultation with the Technical Officer it was decided that the 

model was more reliable if it used the old tiles. 

The temperature fields in the model were revised to allow a more conect venial 
temperature structure to be included. This for the first time has permitted the model to 

calculate accurate storm-time temperatures. 

Extensive comparisons were made between the storm versions of the model and data 
from the DE and AE satellites. The results are presented in the Technical report. Because 
of these comparisons corrections have been made to the VSH model. These changes are 
included in the delivered model. 
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Task Report on tasks 4, 5, and 6. 


Nine tasks were listed in the three year period of this grant. It is the second three ol 
these tasks that are of relevance here. These tasks are: 


4) Further checks will made to assess the importance of the Y component of the 
Interplanetary ; Magnetic Field on the neutral winds, temperatures, composition and 
densities. A final decision will be made as to whether the gains made in terms of accuracy 
are more important than the losses due to reduced efficiency. 


5) The importance of variations of the Y component of the IMF during storms will be 
assessed. Changes in the Y component during these times are likely to be more important 
than changes during quiet geomagnetic times , so this separate study is necessary. 


6) If needed, the effects of the Y-component if the IMF on the neutral thermospheu 
during geomagnetic storms will be included in the VSH model. 


During the Pi's visit to Marshall in November the Technical Officer (Jett Anderson i 
stated during discussions that the work relating to IMF effects in the VSH model were not 
of great interest to the Scientific and Technical Staff at Marshall because their main 
problems involved the prediction of thermospheric effects on satellite's in low inclination 
orbit. Instead the Technical Officer (Jeff Anderson) stated that the Staff at Marshall 
would be interested in a system of quick checks to be able to see what the historical 
densities were at a particular site that a satellite might travel through. He indicated that 
the PI should decrease the scope of tasks 4 through 6 without indicating how this was to 
be done. In addition the limitation of funds (only the first half of the second year's funds 
were ever allocated) prevented the employment ot more staff to carry out the woik. 
Despite this w'ork was completed on tasks 4 through 6. It was found that the IMF etteci- 



on composition and neutral densities were not of sufficient importance (being restricted to 
high latitudes) to the mission of Marshall to be included in the model. The main reasons 
for not including these effects are: they are relatively unimportant (see the technical 
report); and that a far greater computer overhead would result from their inclusion. 



PREFACE 


This manual is divided into four sections. Part I provides an overview ot the capabilities ot \ SH 
This part is oriented towards all users of VSH and provides the minimum information needed by an end 
user. Part II details how a FORTRAN program would be written to access the VSH routines. It include^ 
examples of main programs. Part III is a technical description of the underlying atmospheric physics and 
mathematical equations used in the VSH model. It is most useful for space scientists and others who are 
interested in the underlying science that went into the development of the VSH model. The tinal section 
directed towards a systems analyst or systems administrator. It provides information on storage, run-time 
requirements, and installation. 



